wd<- "/Users/chengg/Google Drive/EPICON/Mycobiome/Fungal ITS/statistic/Total.fungi"
setwd(wd)

rm(list = ls())
load("EPICON.data.preparation.RC.bNTI.ted.2019.04.19.Rdata")
library(vegan)
library(ggplot2)
library(colorRamps)
library(ape)
library(umap)
library(plyr)
library(plotly)


set.seed(315)
fung.pc<-pcoa(vegdist(decostand(fung.rar,"hellinger"),"bray"))
envpc<-data.frame(fung.pc$vectors, env)

pc1=round(fung.pc$values$Relative_eig[1], 3)*100
pc2=round(fung.pc$values$Relative_eig[2], 3)*100

envpc$Habitat<-factor(envpc$Habitat, levels = c("Leaf", "Root", "Rhizosphere", "Soil") )
envpc$Treatment1<-factor(envpc$Treatment1, levels = c("Control", "Pre_flowering_drought", "Post_flowering_drought"), labels = c("CON", "PRE", "POST") )

df <- envpc[,c("Axis.1", "Axis.2", "Habitat")]
find_hull <- function(envpc) envpc[chull(envpc$Axis.1, envpc$Axis.2), ]
hulls <- ddply(df, "Habitat", find_hull)

p1<-ggplot(data=envpc, aes(x= Axis.1 , y= Axis.2,group=Habitat))+ 
  geom_point(aes(shape=Habitat, colour=factor(TP)),size=2) + 
  geom_polygon(data=hulls, alpha=.1)+
  labs(x = sprintf("PCo1 (%.1f%%)", pc1), y = sprintf("PCo2 (%.1f%%)", pc2))+
  scale_shape_manual(values=c(0, 6, 3, 1))+
  scale_colour_manual(values=blue2red(18))+theme_bw()+
  guides(color=guide_legend(title= "Week"),
         shape=guide_legend(title= "Compartment"))+
  theme(legend.title = element_text(colour="black", size=15, face="bold"),
        legend.text = element_text(colour="black", size=15, face="bold"),
        axis.text=element_text(size=10,face="bold"),
        axis.title=element_text(size=15,face="bold"))
p1

ggplotly(p1)

library(plotly)

plot_ly(x=envpc$Axis.1, y=envpc$Axis.2, z=envpc$Axis.3, type="scatter3d", mode="markers", color=envpc$Habitat)

envpc.sub<-envpc[envpc$TP<2 &envpc$Treatment1=="CON",]
df <- envpc.sub[,c("Axis.1", "Axis.2", "Timepiont", "Habitat")]
df$TPHabitat<-droplevels(interaction(df$Timepiont, df$Habitat))
find_hull <- function(envpc.sub) envpc.sub[chull(envpc.sub$Axis.1, envpc.sub$Axis.2), ]
hulls <- ddply(df, "TPHabitat", find_hull)


p2<-ggplot(data=df, aes(x= Axis.1 , y= Axis.2,group=Habitat))+ 
  geom_point(aes(shape=Habitat, colour=Timepiont),size=2) + 
  geom_polygon(data=hulls, alpha=.1)+
  labs(x = "PCo1", y = "PCo2")+
  scale_shape_manual(values=c(3,17,16, 15))+
  scale_colour_manual(values=c("blue", "red", "brown"))+theme_bw()+
  guides(color=guide_legend(title= "Week"),
         shape=guide_legend(title= "Compartment"))+
  theme(legend.title = element_text(colour="black", size=15, face="bold"),
        legend.text = element_text(colour="black", size=15, face="bold"),
        axis.text=element_text(size=10,face="bold"),
        axis.title=element_text(size=15,face="bold"))

p2

ggplotly(p2)

subs<-((envpc$TP == 8  & envpc$Treatment1 != "POST") | 
         (envpc$TP == 17 & envpc$Treatment1 != "PRE")) & 
  (envpc$Habitat=="Root" | envpc$Habitat =="Rhizosphere")
p3<-ggplot() + geom_point(data=envpc[subs,],aes(x= Axis.1 , y= Axis.2,shape=factor(TP), colour=Treatment1),size=3,alpha = 0.6) + 
  labs(x = sprintf("PCo1"), y = sprintf("PCo2"))+
  scale_colour_manual(values=c("black", "blue", "red"))+theme_bw()+
  facet_wrap(~Habitat,scales="free", nrow=1)+
  guides(color=guide_legend(title= "Treatment"),
         shape=guide_legend(title= "Week"))+
  theme(strip.text.x = element_text(size = 15,face="bold"),
        legend.title = element_text(colour="black", size=15, face="bold"),
        legend.text = element_text(colour="black", size=15, face="bold"),
        axis.text=element_text(size=10,face="bold"),
        axis.title=element_text(size=15,face="bold"))
ggplotly(p3)

set.seed(315)
fung.umap<-umap(decostand(fung.rar,"hellinger"))
envumap<-data.frame(fung.umap$layout, env)


envumap$Habitat<-factor(envumap$Habitat, levels = c("Leaf", "Root", "Rhizosphere", "Soil") )
envumap$Treatment1<-factor(envumap$Treatment1, levels = c("Control", "Pre_flowering_drought", "Post_flowering_drought"), labels = c("CON", "PRE", "POST") )


p4<-ggplot() + 
  geom_point(data=envumap,aes(x= X1, y= X2, shape=Habitat, colour=factor(TP)),size=2)+
  labs(x = "Umap1", y = "Umap2")+
  scale_shape_manual(values=c(0, 6, 3, 1))+
  scale_colour_manual(values=blue2red(18))+theme_bw()+
  guides(color=guide_legend(title= "Week"),
         shape=guide_legend(title= "Compartment"))+
  theme(legend.title = element_text(colour="black", size=15, face="bold"),
        legend.text = element_text(colour="black", size=15, face="bold"),
        axis.text=element_text(size=10,face="bold"),
        axis.title=element_text(size=15,face="bold"))
p4

ggplotly(p4)